function [err]=fdnoise2s(AN,D1,D2,M)
% fdnoise2s(AN,D1,D2,M)
% 	Parameters:
% 		AN:	number of points to be used in the time seriese shifted by this function
% 			NOTE: this is approximately twice the number that will be generated for the
% 			initial sequence of complex numbers in the Fourier domain.
% 		D1:	First displacement
% 		D2:	Second displacement
% 		M:	length of kernel
% 
% 	Outputs:
% 		err	The fractional error resulting from comparing the power spectrum resulting from
% 			perfoming the shift as a whole to the shift being done in two parts
% 
% 	Purpose of this code:  This code is a modification of the code fdnoise.  It is specifically designed
% 	to perform the same shift performed in fdnoise, but now to do it once in one step and then again
% 	in two distinct steps.  The fractional error that results from comparing the two power spectrum 
% 	generated can be plotted by this code.

hold off;

%% establish variables
D = D1 + D2;

LB = floor(D+M);

SHFTD = floor(D);

if (mod(LB,2) == 1)
    LB = LB + 1;
    SHFTD = SHFTD + 1;
end

N = AN + LB;

x = 1:AN;
H = AN./2;

%% generate noise
a = random('norm',0,1,1,H-1);
b = random('norm',0,1,1,H-1);

c = a+b.*i;
k = fliplr(conj(c));

%% generate sequence to introduce to ifft
s = [0 c 0 k];
s = s.*(sqrt((N-1)./2));

%% IFFT from FD to TD
ts = ifft(s);

%% Plot Power Spectrum of Original Dataset
ps = (c.*conj(c));
xps = 1:(H-1);
%plot(xps,ps,'g');
%hold on;

%% Plot Shifted Time Sequence by D1 + D2 (Interpolation)
tail = zeros(1,SHFTD);
ats = [ts tail];
sts = FDtest(ats,D,M);
ststr = sts(SHFTD+1:end);
stsf = fft(ststr);
pss = (stsf.*conj(stsf))./(sqrt((AN-1)./2));
psss = fftshift(pss);
b = floor((AN+2)./2);
pssst = psss((b+1):AN);
%plot(xps,pssst,'b');
%hold on;

%% Plot Shifted Time Sequence by D1 (Interpolation)
sts1 = FDtest(ats,D1,M);
%sts1 = shiftparts(ats,D1,M);


%% Plot Shifted Time Sequence then by D2 (Interpolation)
sts2 = FDtest(sts1,D2,M);
ststr2 = sts2(SHFTD+1:end);
stsf2 = fft(ststr2);
pss2 = (stsf2.*conj(stsf2))./(sqrt((AN-1)./2));
psss2 = fftshift(pss2);
b2 = floor((AN+2)./2);
pssst2 = psss2((b2+1):AN);
%xps2 = xps + (1./AN);
%plot(xps,pssst2,'r');
%xlim([50,N]);

%% Calculate and Plot Error
err = abs(pssst2 - pssst)./(abs(pssst));
%plot(xps,err,'g');

semilogy(xps,err,'g');
%hold on;
%end
end